Shear yielding of amorphous glassy solids: Effect of temperature and strain rate 
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We study shear yielding and steady state flow of glassy materials with molecular dynamics sim- 
ulations of two standard models: amorphous polymers and bidisperse bennard- Jones glasses. For a 
fixed strain rate, the maximum shear yield stress and the steady state flow stress in simple shear 
both drop linearly with increasing temperature. The dependence on strain rate can be described 
by a either a logarithm or a power-law added to a constant. In marked contrast to predictions of 
traditional thermal activation models, the rate dependence is nearly independent of temperature. 
The relation to more recent models of plastic deformation and glassy rheology is discussed, and the 
dynamics of particles and stress in small regions is examined in light of these findings. 

PACS numbers: 83.60.La, 83.10.Rs, 64.70.Pf 



I. INTRODUCTION 

Deformation processes and plasticity in amorphous 
materials such as metallic or polymeric glasses have re- 
cently received a lot of attention [D, 0, 0, El HI • These ma- 
terials are used in many load-bearing applications. How- 
ever, an understanding of their yield and flow properties 
is hampered by the absence of long-range order and eas- 
ily identifiable mechanisms that mediate the deformation 
such as dislocation motion in crystals. A similar situ- 
ation is encountered in "soft" glassy materials such as 
foams, pastes, and colloidal suspensions, which are also 
characterised by a liquid-like structure and long relax- 
ation times. In fact, it has been suggested recently that 
these very different materials can be viewed as particu- 
lar realizations of a jammed state |(|, which implies that 
their mechanical behavior could be described in a com- 
mon framework. 

Much insight into the mechanical behavior of struc- 
tural glasses has been gained from molecular simulations 
of simple glass-forming liquids or polymers, where parti- 
cles interact through a Lennard-Jones potential. For ex- 
ample, Falk and Langer 3] studied two-dimensional shear 
deformation of a mixture of such particles and found lo- 
calization of plastic events in so-called shear transforma- 
tion zones. Barrat and Berthier studied the steady 
state flow of a similar model and analyzed its relation to 
the fluctuation-dissipation theorem in out-of-equilibrium 
situations. 

In a recent paper , we have studied the onset of shear 
yielding of amorphous polymer glasses under multiaxial 
loading conditions. A pressure-modified von Mises crite- 
rion Q accurately describes the maximum shear yield 
stress as a function of the applied stress for different 
temperatures. However, in these simulations the bulk 



•present address: Laboratoire de Physico-Chimic Thcorique, 
ESPCI, 10 rue Vauq uelin, F- 75231 Paris Ce dex 05, France 
^Electronic address: Joerg.Rottler@jhu.edu 



polymer was deformed at a single constant strain rate. 
The present work investigates the effect of strain rate on 
the stress at the onset of shear yielding as well as on 
the steady state flow stress of glassy materials in simple 
shear. In order to relate the onset of yield to steady shear, 
we not only discuss the macroscopic material response, 
but also perform an analysis of the localized dynamics 
of the stress distribution in the deforming solid. Such a 
program is particulary suited to test predictions of mod- 
els of plasticity 0, S H E E3 and may lead to a deeper 
understanding of the nature of plastic deformation. 

In the following section, we briefly summarize tradi- 
tional and more recent models of viscoplasticity and rhe- 
ology of glassy materials. Sections II I II and IIVI discuss 
the molecular models and the simulation results, respec- 
tively. Section [V] critically reviews how the theoretical 
ideas of Section [Til describe the simulation results. 



II. MODELS OF PLASTICITY AND 
RHEOLOGY 



Rate and temperature dependence: Eyring 
model 



The simplest model that makes a prediction for the 
rate and temperature dependence of shear yielding is the 
rate-state Eyring model 0, of stress-biased thermal 
activation. Structural rearrangement is associated with 
a single energy barrier E that is lowered or increased 
linearly by an applied stress a. This defines transition 
rates of the form 
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where vq is an attempt frequency and V* is a constant 
called the "activation volume." In glasses, the transition 
rates are negligible at zero stress. Thus at finite stress 
one needs to consider only the rate i?+ of transitions in 
the direction aided by stress. The plastic strain rate e p i 
will be proportional to R + , e p i — cR + . Solving for the 
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stress a, one obtains 
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Eq. J2)l contains only a single relaxation time scale and 
predicts an apparent yield stress that varies logarithmi- 
cally with the strain rate and a prefactor that depends 
linearly on temperature. Despite its simplicity exper- 
imental results Q are often fitted to Eq. and the 
value of V* is associated with a typical volume required 
for a molecular shear rearrangement. 



B. Modern approaches 

Modern phenomenological approaches pay tribute to 
the complexity of glassy systems through several exten- 
sions. First, it has been realized that assuming a single 
energy barrier for rearrangements is an oversimplified de- 
scription of glassy materials [l| . One can therefore intro- 
duce a distribution of barriers and add additional time 
scales. Second, any theory that attempts to predict a full 
stress-strain curve must contain some information about 
the internal state of the system as a function of time or 
strain. Extensions therefore consider dynamical internal 
state variables. In the following, we describe particular 
realizations of these ideas. 



1. Shear transformation zone theory 

Falk and Langer used molecular dynamics simulations 
very similar to the present study to identify local plas- 
tic rearrangements under shear. They formulated a the- 
ory of viscoplasticity based on the concept of "shear 
transformation zones" (STZ), bistable (mesoscopic) re- 
gions that transform under shear between ±-states. One 
then considers the dynamics of an ensemble of STZ with 
number density n± on a mean-field level, which deter- 
mines the plastic strain rate 



Ao(R + n + — R-ri-), 



(3) 



where Aq is a constant. A key difference from the Eyring 
model resides in the form of the transition rates R±, 
which are assumed to be free-volume (entropically) ac- 
tivated rather than thermally activated, i.e. 



R± = i?o exp 



v exp [tct/m] 



(4) 



where vq is a characteristic free volume required for a 
STZ flip and /2 a characteristic stress scale required for 
a molecular rearrangement. The role of temperature is 
played by a "free volume" vf per particle. The authors 
motivate this with the observation that in a solid at very 
low temperature, energy barriers should be very large 
compared to thermal energies and thus, as in granular 



systems, thermal activation over these barriers should be 
negligible. The population densities themselves evolve 
according to the rate equation 



h± = R T n T - R±n± + <re p i(A c - A a n±) 



(5) 



where the last term introduces creation and annihilation 
processes of STZ's proportional to the work of plastic 
deformation ae p i. The STZ equations can be solved ana- 
lytically in special steady state situations and otherwise 
solved numerically. They were shown |3| to have both a 
jammed solution for which e p i — and a flowing solution 
once a exceeds a true yield stress a y . The shear rate 
rises linearly as sigma increases above a y as in a Bing- 
ham fluid. Numerical stress-strain curves, hysteresis ex- 
periments and creep tests were shown to be accurately 
reproduced by the model after the adjustment of several 
fit parameters. 

Lemaitre recently extended STZ theory with con- 
cepts from the physics of granular media. He treated the 
free volume parameter Vf as a dynamical state variable 
and proposed the following time evolution: 



-Ri exp 



Vf 



A v ae p i. 



(6) 



This expression is motivated by slow density relaxations 
in granular materials that decrease the free volume. The 
activation factor exp [— Vi/vf] describes the probability 
for volume fluctuations larger than a characteristic vol- 
ume Vi. Note that the "activation barrier" for com- 
paction Vi differs a priori from the barrier for shear trans- 
formation vo- The second term refers to creation of free 
volume (shear induced dilatancy) again due to plastic de- 
formation. A " linearized" version of this theory produces 
a power law relation between shear rate and shear stress, 
c ~ £pi K ~ 1 +1 , where k = Vi/vq. The full theory can 
yield more complicated functional forms with a true yield 
stress 01. 



2. Soft Glassy Rheology model 

The soft glassy rheology (SGR) model of Sollich et al. 
[Toj is an extension of a trap model for glasses originally 
proposed by J. -P. Bouchaud ^3 with stress acting as 
an external drive. It was designed to describe the flow 
behavior of foams, dense emulsions, pastes and slurries. 
However, it is very similar to the previous models and 
should also be relevant to the materials of interest here. 
Small volume elements are assumed to yield with a rate 
roexp[— (E — kl 2 /2)/x], where I is the local strain and 
k an elastic constant. The role of temperature is re- 
placed by a "noise" temperature x, which is assumed 
to describe the effect of structural rearrangements in a 
mean- field spirit. Note that the model describes stress- 
assisted yielding as in the other two models, but stress 
enters quadratically and not linearly. Structural disor- 
der is modeled with an exponential distribution of yield 
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energies E, p{E) — exp(—E/x g )/x g , where x g is usually 
set to one. This introduces a distribution of relaxation 
timescales. After yielding of a volume element, a new 
yield energy is drawn from p{E) and the local strain I 
rises again from according to a macroscopic shear rate 
7. The time evolution of the probability P{l,E\t) de- 
scribing the ensemble of volume elements can be obtained 
from a master equation. 

Like the STZ theory, the predictions arising from the 
SGR equations are richer than the simple Eyring model. 
The exponential distribution of traps induces a dynam- 
ical glass transition, and the system exhibits aging for 
x < 1. Analysis has mainly focused on the steady state 
situation under constant shear rate 7, which is the generic 
experiment used to determine the mechanical properties 
of soft glassy materials. Salient predictions are: a New- 
tonian fluid flow a cx 7 for x > 2 and a power law fluid 
a oc j x ~ 1 for 1 < x < 2 In the glassy phase, a 

scaling of the form a — a y oc •j 1 ~ x is predicted. The 
nature of the noise temperature x (not the true thermo- 
dynamic temperature) and the pre-exponential factor Tq 
(thermal) remain largely unspecified. 

3. Microscopic Approaches 

The descriptions of nonlinear rheology and plasticity 
described so far are appealing because of their simplicity, 
but much of the physics is put in "by hand" . In the most 
recent literature, efforts are made to derive the response 
of driven glassy systems from microscopic considerations. 
Berthier et al. consider a generic driven glassy sys- 
tem by calculating the correlation and response functions 
from the microscopic Langevin equations in a mean-field 
approximation. Based on this approach, they suggest a 
"two-time-scale scenario," in which the slow time scales 
associated with structural relaxations are accelerated by 
the drive, while the fast degrees of freedom (phonons) 
remain at the thermodynamic temperature. This con- 
cept can be extended to introduce an "effective tempera- 
ture" different from the thermodynamic temperature. A 
relevant conclusion for the present work is that in their 
calculations for T > T c 1 3J , the slow relaxation time t a 
decreases with the drive, t a ~ 7~ 2 / 3 . Since the relax- 
ation time determines the viscosity, this leads to power 
law shear thinning, a ~ 7™ with exponent n — 1/3. Be- 
low T c , their numerical results indicate that n is only very 
weakly temperature dependent. 



III. MOLECULAR SIMULATIONS 

We perform three-dimensional molecular dynamics 
(MD) simulations for two model glasses. A bead-spring 
model is used for polymers. Beads of mass m interact via 
a conventional 6-12 Lennard-Jones (LJ) potential. All 
results will be expressed in terms of the characteristic 
length a, energy uq and time tlj = \J ma 2 /uq of this 



potential. Unless otherwise noted, the potential is trun- 
cated at r = r c = 1.5a for computational convenience. 
We construct linear polymer chains by connecting ad- 
jacent beads with the finite extensible nonlinear elastic 
(FENE) bond potential [l4|. This polymer model has 
been used extensively to study polymer melt dynamics 
|15| and was also used in our previous study of yield con- 
ditions 0- The number of beads per chain used below 
is usually 256, but the chain length and entanglement 
effects are not important for the small strains considered 
in Section HVAl 

In addition to the polymer, a binary mixture composed 
of 80% A-particles and 20% B-particles without cova- 
lent bonds is also studied. LJ interaction parameters 
were set to the values employed in previous studies that 
aimed at verifying predictions of mode-coupling theory 
for supercooled liquids UM. and studied aging [TtJ or dy- 
namical heterogeneities [lg during the glass transition. 
A very similar system with slightly different parameters 
was used by Falk [3| to study deformation and plasticity 
in amorphous metals in two dimensions. 

Both the polymer and binary mixture models enter 
an amorphous glassy state without crystallization upon 
cooling. For the polymer model, the glass transition tem- 
perature T g w 0.35 ± 0.05uo/kB, while T g is smaller for 
the binary system. These values are for r c = 1.5a, and 
slightly higher values of T g are obtained with larger r c 
hj. Here, we focus on a temperature range between 
T = O.Oluo/^s an d T = 0.3uo/fcs- Unless noted, the 
temperature is controlled with a Langevin thermostat 
(damping rate 1t£j ). The equations of motion are solved 
using the velocity Verlet algorithm with a timestep of 
dt = 0.0075 r LJ . 

The simulation cell contained 32768 LJ beads. Pre- 
vious studies of the yield behavior Q have shown that 
an increase of the system size beyond this size leads to 
slightly lower values of the shear yield stress, but the 
generic behavior is unchanged. In order to minimize sta- 
tistical fluctuations while varying the strain rate, we use 
the same initial state for all runs at a given temperature. 
The glassy states were prepared by a quench from a fluid 
temperature of T = 1.3uo/fes to a glassy temperature of 
T = 0.3uo/fcs at constant volume over a time interval of 
order IOOOtlj. The density was chosen so that the hy- 
drostatic pressure was zero at this temperature. Lower 
temperatures were then reached by cooling at constant 
pressure. 

In studies of initial yield behavior, periodic boundary 
conditions are applied in all directions to eliminate edge 
effects. The original cell is cubic and at zero hydrostatic 
pressure. Tensile or compressive strains are imposed on 
one or more axes by rescaling the periods Li. We use 
true strain rates ej = L~ 1 dLi/dt between e, = 10~ 6 r^j 1 
and £j = 10 _3 T LJ 1 . Note that since tlj ~ 3ps, the strain 
rates employed are much higher than typical experimen- 
tal strain rates. As discussed below, sound propagation 
is too slow for stress to equilibrate across the system at 
the highest strain rates. However, most simulations are 
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FIG. 1: Stress-strain curves for (a) the polymer glass and 
(b) the 80/20 LJ glass at T = Q.Qlu /k B and three different 
strain rates 10 _3 r T ~ I 1 , lO -4 '^ 1 , and 10 _s r£j L (from highest to 
lowest curves). 

slow enough that loading proceeds nearly quasistatically. 

Steady state flow cannot be investigated by deforming 
the simulation box in the above manner, because some 
box periods would soon decrease to a single molecular di- 
ameter. One approach is to use Lees-Edwards boundary 
conditions |2£j ■ Here, we choose a different route and re- 
place the periodic boundary conditions in the 3-direction 
with two rigid walls composed of 2 layers of an fcc(lll) 
crystal for simulations of simple shear. Otherwise, the 
simulation cell has the same dimensions and size, and 
the wall beads are strongly coupled to the sheared glass 
so that slip at the interface is prohibited. By moving 
one wall parallel to the solid at constant velocity v x , a 
steady state shear profile is imposed, and we can mea- 
sure the shear stress r as a function of the average strain 
7 = tv x /h, where h = 32a is the wall separation and t 
is the elapsed time. In these simulations, the Langevin 
thermostat is only coupled to the irrelevant (y) direction 
perpedicular to the flow (x) and velocity gradient (z) di- 
rections pi). 



IV. RESULTS 
A. Onset of shear 

In ref. we analyzed the onset of shear deforma- 
tion in polymer glasses. Results for a wide range 



of multiaxial loading conditions, temperatures and po- 
tential parameters were consistent with the pressure- 
modified von Mises criterion. In this criterion, the 
driving force for shear is the deviatoric stress Tdov = 

((oi - a 2 f + (cr 2 - cr 3 ) 2 + (<73 - ci) 2 ) 1/2 A where the 
Oi are the three eigenvalues of the stress tensor. Rather 
than having a sudden onset at a finite strain, irrecov- 
erable deformation was observed at arbitrarily small 
strains. The most robust definition of the yield stress was 
T dcv> t ne peak value of the deviatoric stress as a function 
of strain. This quantity increases linearly with pressure 
in agreement with the pressure-modified von Mises crite- 
rion. 

The results in ref. were obtained for a single value 
of the strain rate. Here we focus on the effect of vary- 
ing strain rate on the shear yield stress. The initially 
cubic simulation cell is expanded in one direction at a 
constant strain rate e 2 = L~ x dL z jdt, and volume V is 
conserved by maintaining L x = L y = y/V/L z . Fig. 
shows the deviatoric stress averaged over the entire sim- 
ulation cell as a function of strain at three different strain 
rates. As can be seen, the behavior of the bead spring 
model (a) and the binary LJ glass (b) is qualitatively 
similar. These curves also closely resemble experimental 
stress-strain curves for e. g. polymeric glasses 0- 

For all systems the initial response is nearly elastic, i.e. 
the stress rises almost linearly with strain. In the qua- 
sistatic limit, the initial slope gives the elastic modulus. 
Results for strain rates of 3 x ICpVjJj 1 and below collapse 
onto this quasistatic behavior. As the strain rate rises to 
lO" 3 -^ 1 and beyond, the initial slope grows. The rea- 
son is that stress is no longer able to equilibrate across 
the system. It is well known that the elastic modulus 
of a heterogeneous system is over-estimated by apply- 
ing a uniform strain (the Voight limit Our algo- 
rithm imposes a uniform strain at each step by rescaling 
the cell dimensions and particle coordinates and the re- 
sulting stress will be too high when the strain rate be- 
comes too large. A rough estimate for the characteris- 
tic time for stress equilibration through the system is 
L/c ~ IOtlj, where c is the speed of sound. When 
the strain is stopped suddenly we find an exponential 
stress relaxation with a characteristic time ~ 15tlj. At 
a strain rate of e'j = 10~ 3 t£} , the time to reach a strain 
of 2% is 20tlj and becomes comparable to the above es- 
timates. We conclude from Fig. that the shear rate 
is slow enough to produce a quasistatic deformation for 
e'j < 10 -3 Tj~j , but not at the highest shear rates. 

As the strain increases, even the curves for lower strain 
rates split apart and saturate at different maximum 
heights T^ ev . The maxima are broad and centered at 
strains of about 6% in the binary LJ glass and 8% in the 
polymer. The fine structure on the curves corresponds to 
individual plastic yield events that are discussed further 
below. This structure decreases with increasing system 
size and temperature as the fraction of the system in- 
volved with typical events decreases. Results for larger 
systems were consistent with those shown here, but could 
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FIG. 2: Rate dependence of the maximum shear yield stress 
for (a) the polymer glass and (b) the 80/20 LJ glass. The tem- 
perature decreases from T = 0.3«o/&b (bottom data points) 
to T = 0.01 uo/fes (top data points) with intermediate values 
of T = 0.2uo/fcs and T = 0.1 ito/fts. Also shown are fits to 
a logarithmic rate dependence T? ev = to + s In e (solid) and a 
power law r? = to + re n (dashed) with n = 0.2. Fit values 
of s were (a) 0.028 ± 0.003 and (b) 0.022 ± 0.02. 



not be extended over as wide a range of shear rates and 
other parameters. 

Fig. [2] summarizes values for r^ ev obtained from the 
maximum of curves such as those in Fig. ^ as a function 
of strain rate. Four different temperatures were studied, 
with the highest value, T = 0.3 Mo/fee, close to T g and 
T = 0.01 Mo/fcs far away. This covers a much wider range 
than usually explored in experiments. For a given tem- 
perature and most strain rates, both models give nearly 
straight lines in a semilogarithmic plot. The salient ob- 
servation that can be made in this figure is that all curves 
are nearly parallel, and temperature merely changes the 
offset value on the stress axis. 

As suggested by the Eyring-model, an obvious way to 
describe the rate dependence is through a logarithmic fit 
of the form r% ev = To + sine (solid lines). The obser- 
vation of parallel curves then implies a nearly constant 
prefactor s in front of the logarithm, and the variation 
with temperature is described by To. Often the rate de- 
pendence of the shear yield stress of complex materials 
is also described with a power-law added to a constant, 
i. e. r^ ev = r + re" (dashed lines). As can be seen in 
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FIG. 3: Variation of the shear yield stress with temperature 
for (a) the polymer glass and (b) the 80/20 LJ glass. The 
rates vary between 10 _3 t i ]j L and 10 -6 r£] r 1 as in Fig. [21 (from 
top to bottom). Linear fits to the data points for a given rate 
are also shown. 



Fig. [21 such fits provide an equally good description of 
the data. Due to the small variations of r^ ev , a determi- 
nation of the exponent n via best fits is very unreliable. 
Since the curves for different T are nearly parallel, there 
is no reason to expect n to vary. The curves in Fig. [21 
are for n = 0.2, which gave the smallest variation of the 
prefactor r (< 10%) with temperature, but other choices 
between 0.1 and 0.3 are also possible. Since r > for 
T < 0.2uo/kB, the data is not consistent with a pure 
power law. 

Neither functional form provides a good fit to results at 
the highest strain rate 10 _3 Tlj. We have already shown 
that the elastic behavior changes at this shear rate be- 
cause stress can not equilibrate throughout the system. 
The plastic response will also be affected by the dynam- 
ics of stress distribution, and new behavior is expected 
to set in at this rate. Therefore, these points were not in- 
cluded in the fits shown in Fig. [3 Note that experiments 
are always at much lower strain rates, while most simu- 
lations of shear have been done at e'j = 10~ 5 to 10 _1 r£j • . 
A clear distinction between a logarithm and a power-law 
dependence would thus only be possible by measuring the 
yield stress at lower rates, but unfortunately the simula- 
tion time becomes prohibitively large. 

The fact that the curves of Fig. [21 are parallel implies 
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that replotting the shear yield stress as a function of tem- 
perature also leads to parallel behavior. Fig.Q|shows that 
the trend of t^ cv with T is in fact linear, and the slope 
is nearly rate independent. The rate merely changes the 
offset value. 

In order to investigate whether the rate dependence is 
in any way influenced by the methodology, we have ex- 
panded our studies and considered other model param- 
eters, different thermostat methods and loading states. 
In particular, we obtained data analogous to that shown 
in Fig. El for a longer range of the Lennard- Jones poten- 
tial r c — 2.2a, Nose-Hoover and Langevin thermostats 
with different rates, and uniaxial strain as opposed to 
the volume-conserving shear in Fig. ^ We also consid- 
ered a system that was quenched into the glassy state 10 
times faster than the other systems to investigate cooling 
rate effects. All these different situations show essentially 
the same robust scenario: the rate dependence is nearly 
independent of temperature with a constant offset that 
decreases linearly as T increases. The slope s is gener- 
ally unaffected within the noise, but is roughly twice as 
large for uniaxial strain, and slightly larger for rapidly 
quenched initial states. Again, power-law fits could also 
be used to describe the data, and a fixed exponent suffices 
to describe all temperatures. 

B. Steady shear 

Most rheological measurements of soft glassy materials 
focus on the steady state stress and not on the initial 
transient maximum of the stress-strain curve. Thus it 
is interesting to compare the rate dependence for steady 
state shear to the above results. The steady state can not 
be studied with the entangled bead-spring model used 
above H5 , because the shear stresses at accessible strain 
rates would soon break covalent bonds. However previous 
studies with short chains and the binary system 0, 
13 show similar trends with decreasing T. At high T 
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FIG. 4: Stress strain curves for simple shear of the binary 
LJ solid at T = 0.1 ito/fcs for three different shear rates 7 = 
0.01r£j , 7 = O.OOItj-j 1 and 7 = 0.0001^/ (top to bottom). 



(> 0.7wo/fcs for the binary system), Newtonian behavior 
is observed up to the highest practical shear rates. At 
lower temperatures, where the liquid is in a supercooled 
state, there is a crossover from Newtonian behavior at 
low rates to power law shear thinning a ~ 7™ at high 
rates. The crossover rate goes to zero as T decreases to 
T g , and at lower temperatures the shear stress is nearly 
independent of strain rate [23j| . 

We focus here on the binary (80/20) LJ system at 
T < 0.3uo/fcs, varying T from just above T g to well 
below. As described in Section II I II shear was imposed 
by confining the system between walls separated by h 
along the z-direction and translating one at a fixed speed 
v x . Fig. 01 shows typical behavior of the stress as a func- 
tion of the average shear strain 7. There is an initial 
maximum at 7 ~ 0.1 that corresponds to the transient 
behavior studied in the previous section. The stress then 
decreases rapidly to a steady state plateau value. This 
strain softening becomes more pronounced with decreas- 
ing temperature. 

The average value from the steady state region is plot- 
ted vs. shear rate in Fig. [3] At T = 0.3uo/&b (lowest 
curve), the solid exhibits glassy behavior at high rates, 
but is Newtonian at lower rates, which indicates that the 
temperature is still above T g . For all temperatures be- 
low T g the curves are nearly parallel, as for the results 
for the transient maximum shear stress in Fig. [3 A loga- 
rithmic rate dependence may be fitted at small rates, but 
a power-law added to a constant clearly provides a much 
better fit over the whole range of rates. The fits for the 
temperatures T < 0.2uo/fcs shown in Fig.OOuse n = 0.3, 
for which the prefactor r to the power-law deviates from 
unity by less than 10%. Our results for the glassy state 
are very similar to observations reported by Varnik et 
al. [2^|, who used the same model with a larger cutoff 
value r c = 2.5 a for the LJ potential and only considered 
T = 0.2u /k B - 

The rate dependence of the steady shear stress is very 
similar to that of the yield stress (Fig. [2J. In both 
cases, there is a rapid rise above the logarithmic fits 
at rates of 10~ 3 tlj and above. Although in Fig. [S] the 
system has had time to reach a steady state, one may 
wonder whether there is still a change in behavior near 
7 = 10~ 3 Tlj . At higher rates the stress may not relax 
between the local yield events discussed below, leading 
to a more rapid rise in mean stress. The results were 
not extended to higher shear rates because the tempera- 
ture in unthermostatted directions begins to rise slightly 
above the set temperature even at 7 = 10~ 2 r£j . 

A plot of the shear stress vs. rate as in Fig. [5] assumes 
that the local strain rate is equal to the average implied 
by v x /h. Fig. shows the velocity profile for two differ- 
ent temperatures and several different shear rates. The 
first few layers always move with the walls due to the 
strong coupling. At T = 0.2uo/kB, the profile in the 
center region is homogeneous and here the assumption 
that v s /h equals the shear rate is satisfied. However, at 
T = 0.01 wo/fcs shear is localized in the lower 60% of the 



7 



0.8 
0.6 
0.4 
0.2 





10" 6 10" 5 10" 4 10" 3 



10" 



FIG. 5: Rate dependence of the steady state shear stress for 
the 80/20 LJ glass at 4 temperatures T — 0.3 uo/ks (bottom), 
T = 0.2u o /fcs, T = 0.1it o /fc s ,and T = 0.01u /k B (top). 
Also shown are fits to a logarithmic rate dependence Td cv = 
to + sine (solid) and a power law Tdcv = to + re n (dashed), 
where n = 0.3 for all curves. 



simulation cell, and the upper part of the glass moves 
at the constant wall speed. This is a clear indication of 
shear banding. The local shear rate is now larger than 
the average value, but since the change is only by about 
a factor of two, the data points in Fig.^are not dramat- 
ically affected. 

In their simulations at T = 0.2mo//cb, Varnik et al. 
|26j | also found a transition from homogeneous flow to 
shear banding as the shear rate dropped below 10~ 3 r£} . 
In our simulations, shear is homogeneous at this tem- 
perature, but the shorter cutoff in our model implies a 
slightly lower T g . The increasing amount of shear local- 
ization with decreasing temperature is consistent with a 
stress peak that becomes more pronounced as T is low- 
ered (see Fig. . The region with negative slope on the 
stress-strain curve signals a mechanical instability and 
shear localization. This localization is inhibited in the 
simulations that were used to study the transient stress 
maximum, since they used periodic boundary conditions 
in all directions. Varnik et al. j2(| noted that simu- 
lations with Lees-Edwards boundary conditions [23 also 
suppressed shear banding, but that the shear stresses ob- 
tained from simulations with walls and periodic bound- 
ary conditions are similar. 



Analysis of the dynamics of the local stress 
distribution 



One of the great advantages of molecular simulations is 
the ability to reach beyond a measurement of the macro- 
scopic response function and in addition obtain infor- 
mation about the local dynamics and rearrangements in 
nonequilibrium situations. Such information is essential 
to construct a clear picture of the underlying microscopic 
processes. 
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FIG. 6: Rescaled velocity profiles for 3 different shear rates 
7 = IO-'tu 1 (A), 7 = 10-V-j 1 (*), 7 = IO-V-j 1 (■) at (a) 
T = O.OlMo/fcs and (b) T = 0.2u /k B . The velocities are 
rescaled by the wall velocity v x and the positions by the total 
separation h of the walls. 



Here, we follow this route by decomposing the total 
simulation cell into small volume elements of size 7 — 8a 3 
and measuring the local stress tensor in those regions. 
Since the density is close to la~ 3 , these regions contain 
typically 7-8 particles. A particle number of that size 
should be sufficient to constitute a locally transforming 
region as envisioned in several of the theoretical models 
described in Section [n] Some of the models suggested 
that local rearrangements replace thermal noise as the 
source for activation over barriers, which motivates a 
study of the stress changes experienced by local regions. 

Our results for general stress states have shown that 
the deviatoric shear stress T<j ev is the relevant stress ten- 
sor invariant that describes shear deformation. We there- 
fore calculate the local change Acr-y in the stress tensor 
in every volume element during a small time interval and 
then study the size distribution of changes in the scalar 
variable ATd GV calculated from Acr^- (note that Athov is 
always positive, since it refers to the deviatioric stress of 
Aery) . In general, one expects ATd ov to fluctuate even 
in the undriven case, and one should expect to find a 
stationary distribution of stress jumps over not too long 
timescales. Since the glass is out of equilibrium, it is by 
construction not stationary and will exhibit aging phe- 
nomena, etc. Such long timescales, however, are not ex- 
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FIG. 7: Size distribution of jumps Ardcv for the unstrained 
binary LJ glass at four temperatures T = O.Oltio/fes (solid 
line), T = 0.1u /k B ,T = 0.2u /k B and T = 0.3u /k B (long 
dashed line, T increases from left to right). Ard cv was cal- 
culated for a time difference of 7.5tlj. The distributions are 
stationary over ~ 10 tlj. 
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FIG. 8: Size distribution of jumps in t^ cv for the binary LJ 
glass sheared at 7 = 10 -4 Tj7j for two temperatures T — 
0.01wo/&b and T — 0.3uo/k B - Different curves belong to 
strains of 0%, 2.5%, 4.5% and 7% (see Fig. [TJ. Also shown 
as thick solid lines are the steady state distributions at the 
same rate as well as exponential fits (dashed lines) to the tail 
of those distributions. 



plored here. 

Fig. [7] shows an example of such stationary distribu- 
tions for four different temperatures. The jumps corre- 
spond to an average change in local deviatoric stress in 
small regions over a time difference of 7.5tlj. The nar- 
rowest distribution is found at the lowest temperature 
T = O.OImoAb and the distribution widens as the tem- 
perature increases to T = 0.3uo/A:_b- The cause for the 
stress jumps is obviously thermal motion of the particles, 
and only small excursions about their positions can occur 
in the glass (cage effect). The part of the distribution at 
small values of ATd 0V can be fitted to a Gaussian, but 
deviations become visible at large values. 

We are now able to determine how the distribution of 



stress jumps in the undriven system changes when the 
system is macroscopically strained. Fig. |H1 shows several 
distributions at the highest (T = 0.3uo/A;b) and low- 
est (T = O.Oluo/fcs) temperatures considered. ATd 0V 
was calculated for the same time difference 7.5tlj. For 
each temperature, four distributions are shown, which 
correspond to four different strains between zero and the 
maximum shear yield strain. Also shown for comparison 
are the distributions in steady shear. 

Several observations can be made in this figure. Note 
first that at T — O.Suo/ks, the unstrained distribution is 
not changed dramatically under external driving. Only 
the tail at large jumps gets modified. This temperature 
is above T g and the shear rate is low enough that the 
deviation from Newtonian behavior is small. This is in 
sharp contrast to the situation at T — O.Oluo/ks, where 
strain produces a dramatic increase in the number of 
large jumps. The tail of the distribution can be fit to an 
exponential form exp(— Ardov/Tc) that extends to larger 
stresses as the strain increases to the peak in the stress- 
strain curve. The steady state stress results are close 
to those at a transient strain of 4.5%. The steady state 
probability curves decay with the characteristic stresses 
t c = 0.18u /a 3 (T = O.Oluo/M and r c = 0.12w /a 3 
(T = 0.3u /k B ). 



V. INTERPRETATION AND COMPARISON TO 
MODELS OF PLASTICITY 

The above results for the strain rate dependence of 
shear yielding offer an opportunity to test the theoret- 
ical models described in Section [H] The Eyring model 
Eq. J5J predicts a logarithmic dependence of t^ cv , and 
the prefactor s is given by ksT/V* . Since we found typi- 
cal values of s between 0.02 and 0.03 for all temperatures, 
our result can only be reconciled with the Eyring model if 
one allows for huge variations of the "activation volume" 
V* between 0.3 a 3 and 10a 3 . V* is a phenomenological 
fit parameter, but is typically interpreted as a character- 
istic volume for a local shear event. It should then be 
at least of order the volume per particle, i.e. of order a 3 
or larger. The implied linear change in V* with T would 
also be inconsistent with the observed linear temperature 
dependence of r^ ev . From Eq. J2J) the stress would vary 
as E/V* ~ E/T, which is inconsistent with the data. Be- 
sides, the usefulness of the Eyring expression is reduced 
if V* is allowed to be temperature dependent. 

Modern theories offer a different interpretation of the 
rate dependence. Lemaitre's extension Q of the STZ 
theory [3j to include free volume relaxation produces 
complex rate dependence that could describe our results. 
However, the issue of temperature dependence has not 
been addressed in this approach and we have not exam- 
ined the issue of free volume relaxation here. 

A pure power-law shear thinning behavior was pre- 
dicted in the calculations of Berthier et al. • Starting 
from very different considerations, these authors found 
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that the exponent n should have a generic value of 1/3 
at a transition temperature T c . Below T c , the exponent 
should in principle be temperature dependent, but their 
numerical results indicated only a very weak dependence. 
Both results are consistent with our findings near T g , but 
not at lower temperatures where we find a constant plus 
power-law behavior. 

Below the glass transition temperature, the SGR 
model ^3 predicts such constant plus power-law behav- 
ior of the flow curve with an exponent of the form 1 — x, 
where x is an effective noise temperature in units where 
the glass transition temperature x g = 1. In steady shear, 
we found good agreement with this functional form and 
an exponent n = 0.3 ± 0.1. The authors of the SGR 
model go to great length in providing an interpretation 
for x. They argue that in a sheared state, the energy 
for surmounting a barrier for rearrangement is not only 
provided by the thermal energy, but also by the energy 
released from rearrangements elsewhere in the material. 
This energy diffuses through the material and provides 
an effective thermal bath in a mean field sense. The en- 
ergy released in such a rearrangement must therefore be 
of order of the typical yield energies, which implies that 
x should be of order unity. When the yield energies are 
much larger than typical thermal energies, x will be inde- 
pendent of the true thermodynamic temperature T. Our 
finding of an exponent independent of T could thus be 
rationalized in this framework. 

Our analysis of the distribution of local stress jumps 
could lend additional support to the concepts behind the 
SGR model. Changes in the local deviatoric stress as 
shown in Fig.[S]can activate yield events. We found that 
the distribution of stress jumps could be fitted to an ex- 
ponential distribution with a characteristic decay stress 
that varied less than 30% as T changed from 0.3uo/^b 
to O.Oliio/fc-B- At small values of ATd ev , the distribu- 
tion retains the equilibrium form. This result suggests 
that the drive generates internal dynamics on a common 
scale despite very different thermodynamic temperatures 
and might provide a more microscopic justification of 
the background "effective noise temperature" proposed 
in the SGR model. 

The shear energy released by a yielding local region 
should indeed trigger additional yield events at other lo- 
cations in the material. However, the ensuing dynam- 
ics may be more complicated than suggested by mean- 
field approaches. Many models of material breakdown 
include load redistribution mechanisms. In fiber bundle 
models j2^|, for instance, one finds avalanche behavior 
that precedes total failure. The avalanche size distri- 
bution follows a power law. When shearing a material, 
however, a yielded region will typically reemerge with a 
different yield energy as assumed in the SGR model. In 
ref. p9j , the simpler situation of an elastic network with 
random breaking thresholds of the links was subjected to 
an external drive. Damaged elements were replaced with 
new ones with a different breaking threshold. Power-law 
behavior in the size and duration of failure events was 



found. Recent work on molecular [30j as well as sim- 
plified models [3l| shows that load redistribution in the 
yielding glass mediated by clastic interactions between 
rearranging regions can lead to avalanche behavior. 



VI. CONCLUSIONS 

The (transient) maximum deviatoric shear stresses for 
two model glasses, binary LJ mixtures and polymers, 
were shown to exhibit strikingly similar trends with rate 
and temperature. The rate dependence was remarkably 
constant from ~ T g to temperatures 30 times lower. The 
entire flow curve shifted linearly to lower stress as T in- 
creased. At low rates, the rate dependence of the peak 
stress could be described by a logarithm or constant plus 
power-law, where the prefactor and/or exponent did not 
vary with temperature. A more rapid rise in stress was 
observed when strain rate was increased to 10~ 3 r£y . De- 
viations in the elastic response also set in at this strain 
rate, indicating that stress could not relax throughout 
the system. This is not surprising given that the time to 
strain to 1% is comparable to that of sound propagation 
back and forth across the system. 

The stress for steady shear flow was calculated for the 
LJ mixture. Curves for different temperatures were also 
nearly parallel, shifting rigidly to lower stresses with in- 
creasing temperature. Below strain rates of 10~ 3 tlj, the 
dependence of the flow stress on rate could be described 
by a logarithm. The entire flow curve could be fit by a 
constant plus power law with a temperature independent 
exponent n — 0.3 ± 0.1. However, as for the transient 
stress, the stress begins to rise more rapidly at 10 -3 Tj7j . 
It is interesting to ask whether the emergence of a power- 
law vs. logarithmic behavior is related to the overlap of 
shear rate and stress equilibration timescales. Although 
the system is in steady state, regions of a few atoms un- 
dergo a substantial yield event after of order 1% strain 
increments. At strain rates of lO -3 ^ 1 and above, the 
time between these yield events will be too short for stress 
relaxation in our system. This point should be kept in 
mind when comparing molecular simulations to experi- 
ments at lower strain rates. 

The Eyring model was shown to be incompatible even 
with logarithmic fits to the shear stress at low rates. It 
predicts that the prefactor of the logarithm should scale 
as kT/V* , while the observed prefactor is nearly inde- 
pendent of temperature. Fixing this discrepancy by re- 
quiring V* to scale linearly with T leads to unphysically 
small values of V* and is inconsistent with the linear 
drop in yield stress with temperature at a given rate. 
The Eyring model has been very helpful in analyzing ex- 
perimental data, but typically over a narrow range of 
temperatures. It would be interesting to extend these 
measurements to very low temperatures and to higher 
shear rates. Note that typical room temperature exper- 
imental values for V* in polymers correspond to 3 or 4 
repeat units [8| , which is of the same order as the values 
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we find at T = 0.1 or 0.2u /k B . 

The insensitivity of the rate dependence to temper- 
ature changes indicates that thermal activation is not 
dominant at the rates studied. Analysis of the individual 
particle trajectories offers an explanation for this. Even 
in the limit of zero temperature we find that the trajec- 
tories are exponentially sensitive to changes in rate and 
other parameters. The system does not travel along a sin- 
gle path through the energy landscape at different rates, 
but gets deflected between many possible paths by small 
perturbations. 

Both the STZ and SGR models describe complex dy- 
namics in systems where temperature is unimportant. 
Instead, activation is due to another internal state vari- 
able, either vj or x, that couples to the external drive via 
a feedback mechanism. The original STZ theory gives a 
simple linear rise in stress with rate above the yield stress 
. However recent generalizations Q can produce more 
complex rate dependence like that found here. The SGR 
model predicts a constant plus power-law behavior that 
is also consistent with our simulations and can account 
for a temperature independent exponent n. It remains to 
be seen whether the generating mechanism of structural 
disorder (SGR) or free volume dynamics (STZ) provides 
a more useful description of glassy dynamics. 

An external drive introduces a new timescale into the 
problem that couples to the structural rearrangements. 
It will therefore modify the associated relaxation time 
scales and alter the unperturbed glassy dynamics, e.g to 
stop aging and induce rejuvenation [fj. At zero temper- 
ature, this timescale should be the only one present in 
the system, while other competing timescales will arise 



at finite temperature. Associated with this timescale is 
the concept of an effective temperature 0, that also 
appears in the SGR model. While this picture is clearly 
a useful starting point, the analysis of the jump distribu- 
tion has shown that the internal dynamics of the flowing 
glassy solid is much more intricate and strongly calls for 
extensions of analytical models beyond mean-field con- 
cepts. Such a treatment should describe more accurately 
the load redistribution mechanisms and the propagation 
of released shear energy on a microscopic level. A prelim- 
inary analysis has shown that local plastic events occur 
in avalanches over a wide range of length and time scales. 

Finally, we note that the above analytic models of 
glassy rheology and viscoplasticity are all scalar models 
that cannot describe self-organization on larger length 
scales such as shear bands. Tensorial versions of the STZ 
theory are being discussed to address banding and neck- 
ing |33j| . These are promising approaches that should 
continue to benefit from insight gained from molecular 
simulations, in particular the load redistribution mecha- 
nisms of shear yielding zones. 
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